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Preface 


The  purpose  of  this  study  was  to  evaluate  the  Tersoff 
silicon  potential  for  use  in  Molecular  Dynamic  simulations. 
There  has  been  a  great  deal  of  interest  in  silicon  modeling 
in  recent  years,  but  a  successful,  general  purpose 
temperature  independent  silicon  potential  has  so  far  eluded 
researchers.  Although  the  focus  of  this  work  was  looking 
toward  me  1 1 i ng /annea 1 i ng  studies  to  model  the  reversible 
crystal /amorphous  transition,  a  general  purpose  silicon 
potential  would  be  useful  in  modeling  the  thin-layer  plating 
of  silicon  (for  electrical  component  design)  and  a  host  of 
other  interesting  char ac ter i st i cs . 

While  the  results  of  this  particular  study  were  not 
favorable,  the  general  feeling  among  many  researchers  is 
that  the  Tersoff  style  of  potential  holds  a  great  deal  of 
promise.  Since  a  molecular  dynami c /stat i c  computer  code 
already  exists,  I  hope  this  study  is  continued  by  following 
students.  The  two  appendices  included  at  the  end  of  this 
report  should  help  following  students  to  understand  the 
specifics  of  the  code  begun  by  Paul  Thee  in  his  thesis. 

Three  people  have  provided  invaluable  aid  in  the 
completion  of  this  research.  The  first  is  Capt.  Michael 
Sabochick,  who  was  always  willing  to  set  aside  other  equally 
important  projects  to  explain,  discuss,  or  brain-storm 
various  aspects  of  this  project.  A  second  i ndespensab 1 e 


A 
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source  of  aid  was  Dr.  Jerry  Tersoff,  the  creator  of  the 
potential  evaluated  in  this  study.  He  provided  much-needed 
insight  into  the  potential  function  and  the  results  of 
static  runs  he  had  completed,  and  has  expressed  an  interest 
in  continuing  the  research  dialogue  begun  with  this  study. 
The  third,  and  possibly  most  important  source  of  aid  was  my 
wife  Sue,  who  provided  support  and  understanding  while 
single-handedly  managing  a  very  busy  household.  My  heart¬ 
felt  gratitude  goes  out  to  all  three. 

William  D.  Metzler 
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the  goal  of  this  research  was  to  provide  an  initial 
testing  of  a  new  silicon  potential  proposed  by  Tersoff. 
Molecular  Dynamic  ( MD >  simulations  were  performed  at  four 
temperatures  below  the  melting  point  of  silicon  to  evaluate 
the  new  potential’s  ability  to  reproduce  the  characteristics 
of  diamond-Si.  A  constant  enthalpy  ensemble  composed  of  64 
atoms  was  used  for  the  evaluation. 

Although  the  Tersoff  potential  does  identify  the  diamond 
structure  as  the  lowest  energy  crystal  at  absolute  zero» 
problems  were  noted  at  higher  temperatures.  At  all 
temperatures  above  zero,  the  calculat ional  cell  showed  a 
marked  decrease  in  volume.  The  amount  of  system  shrinkage 
increased  with  increasing  simulation  temperature.  Because 
of  this  shrinkage,  the  Tersoff  potential  is  unusable  in  MD 
s i mu  1  at  ions . 

Several  attempts  were  made  to  isolate  and  identify  the 
cause  of  the  shrinkage,  but  were  unsuccessful.  It  was  noted 
that  as  the  bonding  environment  became  more  favorable 
(stronger  bonding),  the  bonding  radius  of  minimum  energy 
decreased.  Since  the  potential  relies  primarily  on 
coordination  number  to  stabilize  the  diamond  structure, 
it  is  postulated  that  there  exists  an  as-yet  unidentified 
four-bond  configuration  that  provides  a  lower  energy  and  as 
a  result  has  a  shorter  rest  bond  length  than  the  tetrahedral 


diamond  structure., 


EVALUATION  OF  TERSQFF’S  revised  silicon  potential 

IN  A  MOLECULAR  DYNAMIC  SIMULATION  OF  SOLID  SILICON 

^  i_on 
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Since  its  inception  in  1959  (1),  molecular  dynamics 
( MD )  has  become  a  very  popular  research  tool  for  the  study 
of  material  properties  (3  -  5,  7,  11,  14).  Unlike  Brownian 

dynamics  or  Monte-Carlo  techniques,  MD  solves  Newton’s 
equation  of  motion  to  produce  a  strictly  deterministic 
simulation  of  the  trajectories  of  the  individual  atoms  in 
the  model  volume.  Monitoring  the  movement  of  each  atom  as 
it  i-teracts  with  neighboring  atoms  through  a  pre— 
established  interatomic  potential  permits  researchers  to 
s  t  j  c  f  t-e  mechanics  of  effects  observed  in  the  laboratory, 
anp  t"  .  ,ies  information  unavailable  from  tradit  i  ona  1 
-eiSureme-t  me  *  r  c  as  11:"799).  Additionally,  MD  offers 
"ompiere  c:rf:;  of  the  environment,  allowing  the  user  to 
cose-'-e  materia,  behavior  under  conditions  that  would  be 
f ; *  *  i r _ 1 t  or  impossible  to  create  in  a  laboratory  (8:1). 

Before  this  powerful  tool  can  be  used  to  simulate  a 
actual  system,  however,  an  accurate  interatomic  potential 
^unction  for  the  material  must  be  determined.  Past  studies 
ta.e  relied  heavily  on  pair  potentials  (6,  7,  9,  11).  The 


Figure  1:  Diamond  Structure  for  Silicon 


advantage  of  pair  potentials  is  that  they  are  relatively 
simple  and  computationally  less  demanding  than  more  general 
potential  functions)  but  still  allow  accurate  modeling  of 
noble  materials.  However,  one  common  feature  of  all  pair 
potentials  is  spherical  symmetry.  This  makes  them 
unsuitable  for  modeling  materials  that  display  directional 
covalent  bonding,  such  as  silicon.  Pair  potentials  fail  to 
stabilize  the  diamond  crystalline  structure  observed  in 
natural  silicon  (Figure  1),  and  degeneration  to  a  hexagonal 
close  pack  ( HCP  )  structure  is  observed  (3:2001). 

A  significant  amount  of  research  effort  has  been 
expended  in  the  development  of  an  accurate  interatomic 
potential  function  for  silicon.  Several  new  potentials  have 


been  devised,  with  most  being  three-body  potentials  (3,  13, 

15).  These  three-body  potentials  provide  stability  to  the 

diamond  crystal  structure,  but  required  calculations  for 

3 

most  three-body  interactions  scale  as  N  ,  where  N  is  the 

number  of  particles.  The  resulting  computational  demand 

limits  the  practical  size  of  systems  that  can  be  studied. 

The  Biswas  and  Hamann  potential  is  an  exception,  because  the 

2 

separable  form  of  their  function  allows  N  scaling  (3:2001). 

The  Tersoff  silicon  potential  offers  a  notably 

different  approach  to  the  problem  of  silicon  modeling. 

Tersoff  has  avoided  explicit  three-body  terms  in  the 

potential,  resulting  in  a  hybrid  two-body  potential  in  which 

the  attractive  term  is  modified  by  the  bonding  environment. 

Because  the  Tersoff  potential  is  basically  a  two-body 

2 

potential,  the  bond  calculations  scale  as  N  rather  than 

3 

N  .  Since  processor  time  is  a  major  limitation  in  MD 
studies,  this  is  a  very  attractive  feature.  In  addition, 
this  potential  includes  structural  information  that  may 
allow  its  use  in  studies  where  silicon  structures  other  than 
diamond  are  to  be  examined.  The  specifics  of  the  Tersoff 
potential  are  discussed  in  Chapter  II. 

To  date,  Tersoff  has  proposed  two  different  silicon 
potentials,  the  second  representing  a  modification  of  the 
original.  Tersoff 's  initial  silicon  potential  (16)  was 
evaluated  in  a  molecular  dynamic  simulation  by  Thee  (18). 

One  major  problem  with  the  original  potential  identified  in 
Thee’s  work  was  that  the  atom  density  increased  with 


vy 


o 

increasing  temperature  up  to  1000  K,  contrary  to  known 
behavior  of  silicon  (19).  As  a  result,  the  original 
potential  was  deemed  unusable  in  MD  simulations.  Tersoff ’s 
modified  potential  function  (17)  was  investigated  in  this 

work  . 


Statement 

The  purpose  of  this  research  was  to  determine  the 

usefulness  of  the  new  Tersoff  potential  for  MD  simulations. 

This  represents  an  initial  step  toward  the  ultimate  goal  of 

simulating  the  melting  and  annealing  behavior  of  silicon 

undergoing  pulsed  laser  heating.  The  melting/annealing  of 

silicon  is  of  particular  interest  to  the  military  as  a 

survivable  means  of  mass  information  storage.  Because  of 

its  annealing  characteristics,  silicon  has  been  widely  used 

as  an  optical  storage  medium.  Advantages  of  optical  storage 

8  2 

include  very  high  data  density  <>  10u  bits/crn  ),  faster 
access  time,  and  resistance  to  EMP  (electromagnetic  pulse  - 
a  nuclear  blast  effect),  unlike  the  magnetic  storage  media 
now  in  wide  useage.  By  studying  the  mechanics  of  silicon 
annealing,  better  optical  storage  designs  may  be  possible. 


potential,  the  simulation  Mas  limited  to  64  atoms  to 
expedite  the  results,  and  periodic  boundary  conditions  were 
used  to  arrive  at  bulk  properties.  Therefore  edge  effects 
and  boundary  phenomenon  were  not  examined. 


A§§y02ilQD§ 

As  with  any  simulation,  numerous  assumptions  and 
approximations  were  required  to  make  this  problem 
manageable.  Rather  than  being  identified  in  this  section, 
required  simplifications  are  discussed  in  the  relevant 
portions  of  the  following  chapters,  where  the  accompanying 
discussion  should  provide  the  reader  with  a  better 
understanding  of  their  effects. 

Round-off  errors  are  the  single  exception.  The 
majority  of  experimental  runs  completed  in  this  research 
effort  were  performed  on  a  Cray  X/MP  computer.  This 
machine  uses  14  significant  digits  in  single  precision,  more 
than  sufficient  to  keep  round-off  errors  negligible,  even 
over  thousands  of  time-steps.  The  remaining  runs  were 
completed  on  a  VAX  11/780  in  double-precision  mode.  This 
provided  16  significant  digits  in  the  numerical  calculation, 
but  the  increased  processing  time  made  it  impractical  to  use 
the  VAX  computer  for  the  majority  of  the  tests  run. 


Genera_l  Approach 

This  research  effort  was  completed  in  four  stages.  The 
first  stage,  consisting  primarily  of  computer  coding, 
involved  modification  of  Thee’s  existing  molecular  dynamics 
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code  (14)  to  include  the  new  Tersoff  potential  and 


appropriate  derivatives.  Once  coding  was  completed*  several 


test  runs  were  performed  to  insure  the  potential  and 


positional  derivatives  (forces)  were  correct.  This 


validation  of  the  code  comprised  stage  two  of  the  research. 


Stage  three  included  the  full-scale  tests  of  the  silicon 


model  at  temperatures  below  melting  to  evaluate  the  ability 


of  the  simulation  to  reproduce  the  physical  behavior  of 


solid  silicon  in  a  diamond  configuration.  Study  and 


analysis  of  the  changes  in  lattice  structure  as  a  result  of 


heating  concluded  the  research 


Sequence  of  Presentation 


The  final  report  on  this  research  effort  is  divided 


into  six  chapters.  Chapter  I  is  the  introduction,  which 


presents  a  brief  history  of  the  development  of  the  silicon 


potential  and  outlines  the  research  procedure.  Chapter  II 


describes  in  detail  Tersoff ’s  new  potential.  Computational 


techniques  used  to  perform  the  simulation  and  procedures 


used  to  calculate  the  material  properties  are  described  in 


Chapter  III.  Chapter  IV  contains  the  results  of  both  the 


test  runs  used  to  validate  the  code  and  the  evaluation  runs. 


which  were  completed  at  four  selected  temperatures  well 


below  the  melting  point  for  silicon.  Behavioral  analysis  of 


the  potential  is  undertaken  in  Chapter  V.  The  conclusions 


resulting  from  this  work  and  recommendations  for  future 


efforts  are  the  subject  of  Chapter  VI. 
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II.  Tersoff’s  Modified  Silicon  Potential 


Interatomic  Po  t  ent  i_a  l_s  IGeneraU 

In  developing  an  interatomic  potential ,  it  is  generally 
assumed  that  the  potential  energy  at  any  point  in  space  can 
be  expressed  as  a  sum  of  contributions  from  individual  x- 
body  interactions.  For  a  system  containing  N  identical 
atoms,  the  total  potential  energy  U  for  the  system  could  be 
written  (15:5262) 

U  <  1 ,  .  .  .  , N  >  =  2  ui<ri>  +  1  u2<rl’r1> 

i  i  ,  j  J 

i  <  J 

i,j,k  3  I  J  k  “N  1  2  N 

i<  j<k 

where  r  is  the  position  vector  for  atom  i,  and  u  is  the 
x  X 

potential  energy  due  to  x-body  interactions,  which  is  a 
function  of  the  position  of  the  x  atoms.  The  first  term  on 
the  right  hand  side  of  (1)  describes  the  potential  due  to 
the  external  environment,  and  is  of  no  concern  in  this  work. 
Therefore  the  first  applicable  term  in  equation  (1)  is  the 
pair  potential  u^»  which  is  summed  over  all  atom  pairs.  The 
u^  term  describes  the  potential  energy  due  to  all  possible 
three  body  interactions,  and  additional  terms  are  necessary 
for  the  higher  order  interactions.  The  final  term  u^ 
describes  the  energy  due  to  the  single  arrangement  involving 
all  N  atoms  of  the  system. 


incorporating  both  the  two-  and  three-body  terms  of  equation 


iii.  The  basic  weakness  in  these  potentials  is  that  they 
are  generally  fit  specifically  to  the  diamond  structure  of 
silicon,  and  have  very  limited  application  in  simulating 
other  silicon  structures  (5:7363). 

The  Tersoff  potential,  on  the  other  hand,  explicitly 
accounts  for  the  bonding  environment,  relying  on  control  of 
tne  coordination  number  as  described  below  to  stabilize  the 
diamond  structure.  This  approach  offers  the  possibility  of 
a  much  more  widely  applicable  silicon  potential. 


I§r sof  f  2.s  New  SU  icon  Potential, 

The  form  of  the  Tersoff  potential  was  motivated  "by 
intuitive  ideas  about  the  dependence  of  bond  order  upon  the 
local  environment"  (16).  The  general  form  of  the  potential, 
as  previously  mentioned,  is  a  modified  pair  potential  very 
similar  to  the  Morse  potential  (16).  Tersoff ’s  proposed 
two-body  potential  is 

U  =  2  u2(rI1)  =  1/2  2  u2(ri1>  <2> 

i  »  j  J  i  *  j  J 

i<)  j*i 


where  the  pair  potential  U2  is  a  function  of  the  separation 


distance  r..  between  atoms  i  and  j  such  that 

ij 


u 2 * r ij >  =  vru)ta  s”‘D<->rij>  -  Bije,<D<-x2rij,]  <3> 


The  f  (r..)  term  is  the  cut-off  function  used  to  limit  the 

C  ij 


range  of  the  potential,  and  B 


contains  the  bonding 


/* 


t  11  tJ  !.♦  I.tj.l 


Table  I:  Parameter  Values  for  Tersoff’s  Silicon  Potential 


Var i ab  1  e 


Va  1  ue 


Var iab le 


Value 


3264.65  <eV) 


22.9559 


95.3727  (eV) 


3.0  (A) 


1 .40949E-1 1 


0.2  (A) 


4.83810 


3.23940  (/A) 


2.04167 


1.32583  (/A) 


8.80498E-6 


3.23940  (/A) 


environment  information  as  defined  below.  The  variables  A 
and  were  chosen  to  provide  correct  bond  energy  and 
bond  length  for  silicon  diamond.  FCC ,  and  simple  cubic 


structures,  as  well  as  other  static  properties  not 
specifically  identified  by  Tersoff  in  his  communique.  The 


values  reported  by  Tersoff  for  these  variables  (and  others 


to  be  introduced  later)  are  recorded  in  Table  I. 


CyizQff  Function  f  .  The  cut-off  function  provides  a 

-C 


It..  a*, 


2.8  2.9  3. 

Separation  Dista. 


3.2  3.3 


gstroms ) 


Figure  2:  Cut-off  Function  vs.  Separation  Distance 


where  r  is  the  distance  in  angstroms  separating  the  atoms  of 
interest.  The  selected  values  for  R  and  D  are  included  in 
Table  I.  For  a  lattice  constant  of  5.4312  angstroms  (solid 
silicon  at  300  degrees  K),  the  3.2  angstrom  cut-off  radius 
falls  between  the  nearest  neighbor  distance  (2.35  angstroms) 
and  the  next-nearest  neighbor  distance  (3.84  angstroms). 

iD^^r  onmen  t  a.1  Factor  B  .  The  environmental  effect  on 
bond  strength  is  reflected  in  the  B  term.  This  term  was 
the  focus  of  Tersoff’s  modification  of  his  original  silicon 


potential  (16).  Only  the  modified  potential  which  was  the 


Oasis  of  this  study  will  be  discussed  here.  Interested 
readers  are  referred  to  the  original  article  to  compare  the 
two  po  tent i a  1 s . 

The  value  of  Bj i  for  any  bonded  i-j  pair  is  found  by 


ij 

Bu  *  V1  *  b2?j  >'1/C2n) 


(5) 


Values  for  Bq  ,  b,  and  n  can  be  found  in  Table  I.  The  z^j 
term  contains  the  bonding  environment  information,  and  is 
def i ned  by 


"lj 


f  )g(0)expC-CX-  (r.,  -  r.  , 


k*i  ,  j 


ik 


3  ik  ij 


( 6 > 


The  cut-off  function  is  as  previously  defined  and 


where  the  summation  is  over  all  k  atoms  for  the  given  i-j 
pair. 

the  g<0)  function  contains  the  angular  dependence  in  the 
potential  . 


g  ( 0 )  =  1  +•  (  -5  )  - 


2  2 

d  Ccos(0)  h] 


(7) 


Theta  is  the  angle  between  the  i-k  and  i-j  bonds,  and  the 
values  of  the  empirical  parameters  ,  c,  d,  and  h  are 
contained  in  Table  I. 

The  calculation  of  presents  a  possible  problem  in 

computer  simulations.  For  configurations  where  r  >>  r  , 

ik  i  J 

the  value  of  can  become  quite  large  (on  the  order  of 

q 

10  ).  When  this  is  then  raised  to  the  n-th  power  (n  23) 


for  the  computation  of  B 


the  capabilities  of  the  computer 


mav  be  exceeded.  This  can  be  avoided  by  noting  that  when 


2ij  becomes  large,  equation  (5)  can  be  simplified  to 


5ij  .  Bq  b-V(2n)  ,-1/2 


(5a) 


with  negligible  loss  of  accuracy.  This  approximation  was 


used  to  evaluate  B.  .  for  large  values  (>  5)  of  z 

1 J  i  j 


Figure  3  gives  a  visual  representation  of  the  effect  of 


changing  environment  on  bond  strength.  Two  trends  should 


be  noted  in  Figure  3.  First,  the  larger  the  B^j  term  is, 


the  lower  the  bond  energy  (stonger  bonding)  for  any 


separation  distance.  Second,  as  B^j  decreases,  the  most 


favored  (lowest  energy)  bond  length  increases.  Therefore 


both  bond  energy  and  preferred  bond  length  are  functions  of 


the  bonding  environment. 


QQ9Edinatign  Number  vs_.  Bond  Strength.  The  competition 


between  coordination  number  and  bond  strength  is  of  primary 


importance  in  determining  the  bond  environment  as 


represented  by  the  B^j  term.  Increasing  coordination 


requires  more  i-j  bonds,  each  of  which  contributes  to  tht 


overall  system  energy.  Since  bond  energies  are  negative^ 


a  greater  number  of  bonds  equates  to  a  lower  (greater 


negative)  energy  and  therefore  a  more  stable  configuration. 


This  is  the  case  with  ordinary  pair  potentials.  In 


Tersoff’s  potential,  however,  increasing  coordination  number 


also  increases  the  number  of  terms  in  the  summation  of  z  , 

ij 


t't 


F i gure  4 : 


vs .  z 


B 


U 


ij 


which  increases  and  decreases  (see  Figure  4>.  The 

decrease  in  B  results  in  a  weaker  ( less  negative)  bond 
J 

energy  for  each  individual  i  — j  bond*  therby  countering  the 
effect  previously  discussed.  This  attempt  to  model  the 
response  of  s  1  1  1  con  — s i  1  i con  bonding  to  changes  in  the 
environment  provides  the  major  d i f f erent i at i on  between  the 
Tersoff  potential  and  other  potentials. 

^DSyiiC  QfBiDdence.  The  angular  dependence  seems  to 
have  only  a  marginal  effect  on  the  bond  strength.  The  rapid 
increase  in  z  with  increasing  coordination  number  is  the 

^  J 
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Figure  5:  Angular  Dependence  of  g<8>  Function 

primary  tool  used  to  stabilize  the  diamond  structure  within 
Tersoff '5  potential.  The  weak  dependence  of  g<9>  on  the 
bonding  angle  0  (see  Figure  5)  is  a  result  of  Tersoff’s 
empirical  fitting  of  the  potential  to  several  static  crystal 
structuresi  and  offers  the  possibility  of  a  potential 
applicable  to  many  different  silicon  structures. 


III.  Molecular  Dynamic  Model 


The  purpose  of  this  chapter  is  twofold:  1)  to  provide 

the  reader  with  a  brief  overview  of  the  process  of  molecular 
dynamic  simulation  to  allow  better  understanding  of  the 
results  reported;  and  S)  to  explain  briefly  the  calculations 
required  to  convert  the  results  of  a  MD  simulation  into  the 
material  properties  reported  in  the  following  chapters. 

If  a  further  explanation  of  MD  or  statistical  mechanics  is 
desired,  the  reader  is  referred  to  texts  8  and  10  listed  in 
the  Bibliography. 

Qy§D^l§!c!  °f  Molecular  Dynamics  1MD2 

The  purpose  of  the  molecular  dynamic  simulation  is  to 
solve  Newton’s  equation  of  motion  for  each  particle  in  the 
simulation  volume,  providing  a  time-history  of  each 
particle’s  position  and  velocity  over  the  period  of  the 
simulation.  For  a  system  of  N  identical  atoms,  this 
requires  the  solution  for  each  time  step  of  a  set  of  N 
differential  equations  of  the  form 

d  d 

m  —  <  v  )  =  F  -  - - CU  ( r  >  ]  (8) 

1  dt  1  1  dr  N 

1 

where  ,  v^ ,  and  r^  are  the  mass,  velocity,  and  position 
of  atom  i,  respectively,  U ( r^ )  is  the  total  potential  energy 
of  the  system  as  a  function  of  the  N  atomic  positions,  and 
F^  is  the  total  force  on  atom  i. 


Analytical  solution  of  the  set  of  differential 
equations  described  by  (8)  above  would  be  impossible  in  a 
practical  sense.  Instead,  the  Gear  Pr ed ic tor -Cor rec tor 
Algorithm  as  described  by  Haile  (8:17)  was  used  to  provide 
a  numerical  solution.  Because  this  calculation  is  the  heart 
of  the  HD  simulation,  a  brief  explanation  of  the  procedure 
will  be  presented  here. 

Gear  Prfdictgr^Cgrrectgr  Ad  go r ^ t hm  1GPCA2 .  Given  the 


positions  r^  and  first  five  time  derivatives  r^^  ,  r^^» 


r^)  ,  ,  and  r^)  for  each  of  the  N  atoms  at  a  given  time 


t  ,  the  GPCA  determined  the  atom  positions  at  the  next 


incremental  time  t’  =  t  +■  dt  through  a  three  step  process 


of  prediction,  evaluation,  and  correction. 

Pr edic t J.on .  A  fifth  order  Taylor  expansion  was 
used  to  predict  the  position  and  first  five  time  derivatives 


for  each  atom  at  time  t’  based  on  the  values  at  t  .  The 


predicted  value  of  the  n-th  time  derivative  of  the  position 
vector  for  atom  i  at  t’  was  calculated  by 


r  \  5-n  /  dt 

pr  n)cc>  =  I  ;  n+k>(t  > - 

1  k=0  1  0  k! 


(9) 


where  pr 


r<n>(t 


*)  was  the  predicted  derivative  value 


dt 


was  the  time  step  size,  and  r^n)<t  >  was  the  n-th  time 


derivative  of  the  position  vector  for  atom  i  at  t  . 

0 


•  Once  the  position  of  each  atom  was 


predicted  for  t’,  the  resultant  potential  energy  and  the 


force  on  each  atom  was  calculated  using  Tersoff’s  potential 
function  (see  Appendix  A  for  derivation  of  forces).  The 


calculated  atomic  forces  were  divided  by  the  atomic  mass  to 

“(2) 

determine  the  actual  acceleration  r^  (t’)  of  each  atom. 

The  evaluation  of  forces  was  the  most  computationally 
demanding  section  of  the  MD  simulation,  and  several  time¬ 
saving  techniques  were  included  in  the  model  to  minimize  the 
number  of  required  calculations.  The  computer  code  used  in 
this  work  was  a  modification  of  the  code  used  by  Thee,  whose 
thesis  provides  a  detailed  explanation  of  these  techniques 
(18).  The  explanation  will  not  be  repeated  in  this  work. 

a  result  of  the  evaluation  step, 

>  the  actual  acceleration  of  each  atom  at  t’  was  known.  The 

9 

difference  between  the  predicted  and  actual  accelerations 

(2) 

far  each  atom  defined  N  error  terms  dr^  Mt’)  by 

dr^2^(t’)  =  r^2)<t’>  -  pr^2 ^  ( t ’  )  (10) 

These  error  terms  were  then  utilized  to  correct  the 
predicted  positions  and  remaining  four  time  derivatives  of 
each  atom  according  to  the  relationship 

r (n)  ( t  ’  )  =  pr(n)<t’>  +  a  dr^n^(t’)  (11) 

i  i  n  i 

The  a^  coefficients  in  equation  (11)  provide  numerical 
stability  to  the  solutions  of  the  differential  equations  and 
vary  with  the  order  of  the  predictor.  Coefficients  for  the 
V  fifth  order  Taylor  expansion  used  in  this  study  are  listed 
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V 

',vj  Table  II:  Values  of  Coefficients  for  5th  Order  GPCA 


Coef f ic ient 

Value 

Coefficient 

Value 

ao 

3/16 

al 

351/360 

a2 

1 

a3 

11/18 

a4 

1  /6 

a5 

1/60 

in  Table  II.  The  corrected  atomic  positions  and  derivatives 
as  calculated  from  equation  (11)  then  became  the  initial 
values  for  the  next  time  increment,  and  were  stored  to 
provide  the  time-history  needed  for  property  calculations. 

Boundary  Conditions  and  Selection  of  N .  The  selection 
of  the  size  of  the  calculat ional  cell  (number  of  atoms  N) 
used  in  the  MD  simulation  represented  a  compromise  between 
the  greater  realism  of  a  system  containing  a  large  number  of 
atoms  and  the  lighter  computational  requirements  (faster  run 
times)  of  a  smaller  system.  With  small  systems*  however, 
only  a  fraction  of  the  atoms  experience  the  normal  4-atom 
coordination.  As  a  result,  edge  effects  would  have 
overwhelmed  the  bulk  properties  of  interest  without 
compensation  for  the  under-coordinated  atoms. 

Periodic  boundary  conditions  were  used  to  eliminate 
these  edge  effects,  as  is  common  in  MD  simulations  of  this 
type.  By  surrounding  the  ca 1 cu 1  a t i ona 1  cell  with  exact 
images  of  itself,  the  proper  coordination  was  provided 
for  each  atom  in  the  cell.  Figure  6  shows  a  two-dimensional 


SO 
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Figure  6:  Two  Dimensional  Periodic  Boundaries  and  Image 

Interaction  (arrows  indicate  velocities)  (18) 

example  of  periodic  boundaries,  indicating  the  interaction 
of  an  atom  in  the  ca 1 cu 1  a t i ona 1  cell  with  one  image  of 
another  atom.  The  presence  of  these  "images"  set  a  lower 
limit  on  the  size  of  the  ca 1 cu 1  at i ona 1  cell.  To  provide 
stability,  the  c a  1 cu 1  a t i ona 1  cell  had  to  be  large  enough  to 
prevent  any  atom  from  simultaneously  interacting  with  more 
than  one  "image"  of  any  neighboring  atom.  A  more  complete 
explanation  of  the  implementation  of  periodic  boundary 
conditions  is  available  in  Thee’s  previous  work  (18),  and 


interested  readers  are  directed  there  for  more  detail.  For 
the  purpose  of  this  research!  64  atoms  in  a  2x2x2  array  of 
diamond  lattices  was  selected  as  an  optimum  size  for  the 
calculat ional  cell.  This  kept  run  times  acceptable!  yet  was 
large  enough  to  remain  stable. 

RyDDIDS  Id?  Simulation.  The  actual  operation  of  the  MD 
simulation  was  divided  into  four  time  phases. 

1.  Initialization 

2.  Seal ing  Steps 

3.  Resting  Steps 

4.  Property  Calculation  Steps 
Each  time  phase  will  be  discussed  separately. 

iDiiiiiizat ion .  To  begin  a  runi  the  64  atoms 
comprising  the  computational  cell  were  assigned  positions 
corresponding  to  the  lattice  sites  in  a  diamond  crystal! 
and  three  random  velocities  (one  for  each  directional 
component)  were  generated  for  each  particle.  The  random 
velocities  were  scaled  so  that  the  linear  momentum  of  the 
system  was  zero  to  prevent  net  translational  motion  from 
causing  erroneous  temperature  calculations.  The  computer 
code  used  in  this  work  also  allowed  the  user  to  save  the 
positions  and  velocities  from  one  run  and  use  them  as  the 
initial  conditions  for  a  following  run!  minimizing  the 
initial  shock  to  the  system  when  runs  at  high  energy 
(temperature)  are  performed.  This  option  was  seldom  used i 
and  the  results  from  runs  begun  in  this  manner  will  be  so 
annotated.  After  assigning  random  velocities!  the  remaining 
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time  derivatives  were  initially  set  to  zero.  This  caused 
errors  in  the  initial  steps  of  the  GPCA,  but  the  derivatives 
qu .  :  k 1 y  corrected  themselves  as  the  simulation  progressed. 

Scaling  Stegs.  The  purpose  of  the  scaling  steps 
was  to  establish  the  desired  energy  level  in  the  system. 
Since  temperature  is  related  to  the  kinetic  energy  of  the 
constituent  atoms,  the  simplest  method  of  controlling  the 
energy  was  to  scale  the  velocities  of  the  individual  atoms 
so  that  the  desired  temperature  was  achieved.  The 
progression  of  the  simulation  then  partitioned  the  resulting 
total  energy  between  potential  and  kinetic  components.  The 
number  of  scaling  steps  required  to  stabilize  the  system 
depended  on  both  the  amplitude  and  the  level  of  accuracy 
desired  in  the  final  temperature,  but  normally  a  few  tens  to 
a  few  hundred  steps  were  sufficient. 

8§§iing  Stegs.  Abrupt  termination  of  scaling 
subjects  the  system  to  a  mild  shock  ( 8 : SB > .  The  purpose  of 
the  resting  steps  (generally  a  few  tens  to  a  few  hundred) 
was  to  allow  this  shock  to  dissipate  before  property 
calculations  were  begun. 

Prgge^ ty  Calculation  Stegs.  The  property 
calculation  steps  generated  the  actual  data  from  the  MD 
simulation.  The  instantaneous  property  values  over  several 
thousand  time  steps  provided  the  data  base  necessary  to 
calculate  the  bulk  properties  of  the  silicon  model.  The 
actual  calculations  required  to  transform  the  instantaneous 
values  into  bulk  properties  will  be  presented  later. 


EnsemD_l.es  Used  and  lime  St  eg  Size  dt.  In  describing  a 
system  of  atoms,  the  six  state  variables  of  density, 
pressure,  volume,  temperature,  enthalpy,  and  internal  energy 
can  all  be  calculated  if  any  three  are  known.  MD 
simulations  are  generally  identified  by  the  three  fixed 
state  variables.  In  this  study,  two  types  of  ensembles  were 
used,  and  each  will  be  described  separately.  The  choice  of 
ensembles  was  found  to  effect  the  size  of  time  step  that 
could  be  used  in  a  simulation.  Ideally,  one  would  like  to 
use  as  large  a  time  step  as  possible  to  simulate  the  longest 
period  of  real  time  in  the  shortest  number  of  time  steps. 
However ,  at  some  maximum  allowable  time  step  size,  the 
numerical  solution  to  the  differential  equations  become 
unstable  and  the  model  fails.  Time  step  size  and  acceptance 
criterion  for  each  ensemble  is  discussed  below. 

Ensemble .  The  traditional  use  of  MD  was 
limited  to  a  "micro-canonical"  ensemble  where  atom  number, 
volume,  and  total  energy  were  held  constant.  During  the 
course  of  the  simulation,  kinetic  energy  (or  instantaneous 
temperature)  and  potential  energy  should  both  fluctuate 
around  average  values  out  of  phase  with  one  another,  so  that 
the  total  energy  remains  constant.  The  acceptable  amount  of 

variation  in  the  total  energy  was  selected  to  be  1  part  in 

15 

10  over  3000  time  steps.  To  achieve  this  consistency,  1  x 

-  3 

10  J  picoseconds  was  determined  to  be  the  maximum  allowable 
time  step  for  temperatures  up  to  1200°  K,  and  this  step  size 
was  used  in  most  of  the  (N,V,E>  ensemble  runs. 
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The  primary  use  of  the  (N,V,E)  ensemble  was  to  verify 
that  the  Tersoff  potential  and  positional  derivatives 
(forces)  were  correctly  coded  in  the  model.  As  a  secondary 
tool,  the  <  N » V  » E  )  simulation  provided  a  cross-check  of 
results  observed  in  the  (N,P,H>  ensemble  described  below. 

The  results  of  these  cross-checks  are  discussed  in  Chapter 
V  . 

S05§!!?52§*  Since  most  laboratory 
measurements  are  conducted  in  a  constant  pressure/constant 
enthalpy  environment,  it  is  difficult  to  directly  compare 
the  results  of  traditional  <N,V,E)  MD  simulations  with 
effects  observed  in  the  laboratory.  Fortunately,  a 
simulation  method  proposed  by  Anderson  (2)  and  later 
modified  by  Parrinelio  and  Rahman  (12)  allows  simulation  of 
an  (N,P,H)  ensemble.  For  a  detailed  explanation  of  the 
development  of  this  method,  the  reader  is  refered  to  the 
cited  articles.  The  Parrinelio  and  Rahman,  or  P&R,  method 
was  used  extensively  in  this  work,  and  a  brief  discussion 
f o 1  lows . 

To  allow  the  c a  1 cu 1  a t i ona 1  cell  to  change  shape,  the 
cell  boundaries  were  defined  by  a  set  of  mutable  basis 
vectors  a,  b,  and  c,  which  made  up  the  columns  of  a 
transformation  matrix  H  (Figure  7).  The  atom  positions  were 
then  described  by  a  set  of  reduced  position  vectors  s^  so 


that 


Figure  7:  Basis  Vectors  of  Computational  Cell 
for  (N,P,H)  Simulation  (14:62) 


The  components  of  the  reduced  position  vectors 
0  to  1 .  The  volume  Q  of  the  computational  cell 
the  determinant  of  the  H  matrix. 


range  from 
simply 


«  =  || H  f  (13) 

Changes  in  the  structure  of  the  computational  cell  are  then 
accomplished  by  altering  the  elements  of  the  H  matrix  in 
response  to  the  difference  between  the  calculated  internal 
pressure  and  the  constant  external  pressure.  All  runs  for 


this  work  were  completed  with  no  applied  external  pressure. 


The  Hamiltonian  for  this  system  as  derived  in  the  P&R 
article  is  merely  the  enthalpy  of  the  system  plus  a  term  to 
include  the  kinetic  energy  of  the  walls  of  the  computational 
cell.  As  pointed  out  in  the  reference  article,  the  wall 
energy  term  is  negligible,  so  that  the  enthalpy  is  a 
conserved  quantity  in  this  ensemble. 

The  P&.R  article  describes  the  transformations  necessary 
to  solve  Newton’s  equations  in  the  reduced  space,  and  to 
repeat  the  explanation  here  would  provide  no  new  insight. 

The  interested  reader  is  referred  to  the  article  listed  in 
the  Bibliography. 

The  desired  consistency  in  the  constant  enthalpy 

4 

simulation  was  set  at  1  part  in  10  over  3000  time  steps. 
This  was  nearly  achieved  with  a  time  step  size  of  2  x  10”^ 
picoseconds  at  1200°  K,  and  this  time  step  size  was  used  for 
the  majority  of  the  (N,P,H)  ensemble  runs. 

C§l£yI§ii9D  9f  Proper t_i.es 

Many  physical  parameters,  such  as  temperature, 
pressure,  and  internal  energy,  are  dependent  only  on  atomic 
positions  and  velocities.  Therefore  each  possible  atomic 
configuration  (as  defined  by  a  unique  set  of  N  position 
vectors  and  N  momentum  vectors)  has  an  associated  set  of 
these  properties.  Statistical  mechanics  calculates  overall 
ensemble  properties  from  expectation  values  based  on  the 
probability  associated  with  each  atomic  configuration  of  the 


ensemb 1 e . 


However,  the  probability  distribution  function 
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for  each  of  the  infinite  number  of  possible  configurations 
must  be  known  to  perform  these  calculations. 

Fortunately)  the  ergodic  hypothesis  equates  the 
ensemble  "average"  properties  (expectation  values)  to  the 
time  average  of  the  ensemble  properties  taken  over  many 
time  steps  (ideally  #  of  steps  approaches  «,  but  in  practice 
several  thousand  steps  is  sufficient).  Therefore  if  a 
sufficiently  large  data  base  exists,  ensemble  properties  can 
be  calculated  from  a  time-history  of  the  atomic  positions 
and  velocities,  without  knowing  anything  about  the 
associated  configurations.  The  MD  simulation  provides  the 
needed  time-history  data  to  perform  these  calculations. 

The  i nstantaneous  values  of  the  temperature  T,  pressure 
P,  internal  energy  U  and  total  energy  E  can  be  calculated 
from  the  atomic  positions  and  velocities  as  follows:  (8,14) 
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T  =  —  I  v  -  v 
3Nk  l  ±  x 


(  14  ) 
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(  16) 
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In  equations  (14  -  17),  k  is  Boltzmann’s  constant,  V  is 
volume,  and  m  is  the  mass  of  the  particles.  All  other 
variables  are  as  previously  defined. 

Using  the  instantaneous  properties  calculated  from  the 
above  equations,  the  ergotic  hypothesis  allows  calculation 
of  the  ensemble  "average1’  over  L  time  steps  by 

1  L 

<Property>  =  -  I  Proper ty( t  )  (10 

L  n=l  n 

where  Propertylt^)  is  the  instantaneous  value  of  the 
property  at  the  end  of  the  n-th  time  step. 

The  summation  term  in  equation  (15)  is  called  the 
vinal  .  This  is  a  term  that  appears  often  in  the 
statistical  mechanic  calculations.  Because  of  it’s 
s i gn i f i canc e ,  Appendix  B  provides  the  derivation  of  the 
ca 1 cu 1  at i ona 1  form  used  in  this  computer  simulation. 


IV.  Results  of  MD  Runs 


The  MD  runs  completed  in  the  course  of  this  research 
can  be  separated  into  test  runs,  which  were  used  to  validate 
the  coding  of  the  potential  and  its  positional  derivatives, 
and  primary  runs,  conducted  to  evaluate  the  performance  of 
the  potential.  The  results  of  both  types  of  runs  are 
presented  in  this  chapter. 

B§§ylt§  of  Runs 

The  first  set  of  tests  performed  were  static 

o 

evaluations  (a  single  time  step)  conducted  at  0  K  to 
reproduce  the  bond  energies  for  the  rest  bond  lengths 
reported  by  Tersoff  (17)  for  the  dimer,  diamond,  simple 
cubic,  and  RCC  structures.  This  was  accomplished  by  setting 
the  lattice  constant  in  the  model  code  so  that  the  nearest 
neighbor  distance  matched  the  bond  length  associated  with 
the  structure  under  study.  The  number  of  atoms  in  the 
ca 1 cu 1  at lona 1  cell  was  then  adjusted  to  prevent  undesired 
interactions  with  the  "images"  created  by  the  periodic 
boundary  conditions.  The  system  energy  was  calculated  by 
the  MD  program,  and  the  value  was  divided  by  the  number  of 
atoms  in  the  system  to  arrive  at  the  energy  per  atom.  Test 
run  results  are  compared  to  Tersoff ’s  values  in  Table  III. 
The  variance  noted  in  the  energy  for  the  diamond  structure 
is  believed  to  be  a  result  of  a  slight  disagreement  (about 
.0018  angstroms)  in  the  rest  bond  length  for  that  structure. 


Table  III.  Calculated  Bond  Energy 
Dimer,  Diamond,  Simple 

and  Rest  Bond  Lengths  for 
Cubic,  and  FCC  Structures 

Bond 

Energy 

Structure 

Rest  Bond  Length 

( eV/atom) 

<  Angstroms ) 

Current 

Work 

Tersoff 

D  i  mer 

B. 3 132 

-1.3117 

-1.31 17 

D i amond 

E. 3518* 

-4.6305 

-4.6298 

Simple  Cube 

a. 5007 

-4.2871 

-4.2871 

FCC 

2 . 7298 

-4.0818 

-4.0818 

*  it  is  be  1 i eved 

Tersoff’s  used  a  1 

slightly  smaller  value 

The  bond  lengths  in  Table  III  were  found  by  calculating  the 
value  of  B^j  for  each  structure*  then  differentiating  the 
pair  potential  to  determine  the  bond  length  of  minimum 
energy.  Moreover,  all  other  values  match  exactly, 
confirming  the  correct  installation  of  Tersoff’s  potential 
in  the  model  code. 

Static  calculations  could  not  be  used  to  validate  the 
force  calculations  in  the  model,  so  a  single  (N,V,E)  run  at 
300°  K  with  a  time  step  of  1  x  10-^  picoseconds  was 
performed.  A  very  small  system  (4  atoms)  was  used  to 
decrease  the  chance  of  fortuitous  cancellation  of  errors, 
and  the  simulation  was  performed  for  5000  time  steps.  The 
temperature,  internal  energy,  and  total  energy  are  shown  in 
Figure  B.  As  predicted,  the  temperature  and  internal  energy 


are  seen  to  oscillate  about  a  mean  value  out  of  phase  with 


Energy  i  (MJ/ktO 


Figure  8:  Temperature,  Internal,  and  Total  Energy  vs.  Time 
(Instantaneous  Values  for  4  Atom  N,V,E  Ensemble 
at  300  degrees  K,  dt  =  1  femtosecond) 


each  other,  and  the  total  energy  remains  constant.  This 
provides  very  convincing  proof  of  the  correct  evaluation  of 
the  atomic  forces,  because  only  an  exact  balance  of  the 
internal  and  kinetic  energies  will  provide  the  consistency 
observed . 


Primary  Si.mulati_qn  Runs 

Four  primary  runs  were  conducted  using  the  constant 
enthalpy  <N,P,H)  ensemble  to  evaluate  the  applicability  of 
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Table  IV:  Identification  Number /Character ist ics 

of  Primary  MD  Runs 


ID  n 

Des i r  ed 

T  emp 

(°  K  ) 

F  i  na  1 

T  emp 

(°  to 

Scaling 

Steps 

Rest i ng 
Steps 

Proper 

Steps 

l 

300 

303.3 

200 

500 

2500 

s 

600 

612.7 

200 

50 

3000 

3 

900 

853.2 

50 

150 

2800 

4 

isoo 

1  135.8 

50 

150 

2800 

the  Tersoff  potential  in  a  molecular  dynamic  simulation. 

All  four  runs  were  conducted  at  temperatures  well  below  the 

known  melting  temperature  of  silicon  in  a  diamond  structure 

(1685°  K),  so  that  only  solid  silicon  was  modeled. 

To  facilitate  the  discussion  of  the  results*  each  run 

will  be  assigned  an  identification  number.  The  identifying 

number  and  characteristics  of  each  run  are  recorded  in  Table 

IV.  All  the  runs  were  completed  using  a  time  step  of  2  x 
_4 

10  picoseconds.  The  differences  in  the  number  of  scaling 
and  resting  steps  are  a  result  of  attempted  optimization  and 
do  not  effect  the  final  results. 

In  general*  the  results  of  the  evaluation  runs  show  the 
Tersoff  potential  to  be  a  well-behaved  function.  The  state 
variables  of  pressure*  temperature*  and  internal  energy  show 
random  fluctuations  about  their  respective  mean  values  over 
the  period  of  simulation,  as  demonstrated  by  Figures  9  -  11. 
The  enthalpy  values  remain  relatively  constant,  as  expected, 
and  the  instantaneous  values  are  shown  in  Figure  IS.  Only 
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Temperature  (°K) 


Figure  10: 


Instantaneous  Temperature  vs 
for  Primary  Runs 


Energy  ■  (MJ/kg) 


Simulation  Time  ( 


Enthalpy  (MJ /kg) 


Figure  15:  Instantaneous  Enthalpy  vs.  Time 

for  Primary  Runs 
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Figure  13:  Instantaneous  Volume  vs.  Time 
for  Primary  Runs 


Figure  14:  Final  Average  Volume  vs.  Final  Temperature 
of  Primary  Runs 


the  volume  behaves  in  an  unacceptable  manner.  From 
experimental  work,  it  is  known  that  crystalline  silicon 
undergoes  expansion  when  heated,  until  melting  occurs  (19). 

A  quick  glance  at  Figure  13  makes  it  clear  that  this  model 
produces  unrealistic  behavior  as  the  temperature  is 
increased,  showing  a  marked  decrease  in  the  system  volume  as 
the  model  is  heated.  To  emphasize  the  temperature 
dependence  of  this  behaviour.  Figure  14  shows  the  average 
volume  (after  .56  psec  simulation  time)  vs.  the  final 
temperature  with  a  spline  fit  connecting  the  data  points. 

The  volume  is  seen  to  drop  with  increasing  temperature,  and 
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V.  °f  Heat rln^uced  Shrinkage 

The  system  shrinkage  noted  in  the  previous  chapter  is 
even  more  severe  than  that  reported  by  Thee  on  investigation 
of  Tersoff’s  earlier  silicon  potential  (IS).  Thee  reported 
a  similar  decrease  in  system  volume  up  to  about  1000°  K,  at 
which  point  the  trend  reversed  itself  and  the  system  began 
to  expand.  In  this  work,  the  volume  loss  showed  no  such 
reversal  over  the  temperature  range  examined  (300  -  1200 
degrees).  The  purpose  of  this  chapter  is  to  report  the 
results  of  tests  run  to  verify  the  shrinkage  problem  and  to 
attempt  to  determine  the  cause  of  the  shrinkage. 

QonfiLGliif ion  of  Problem 

This  research  effort  utilized  a  molecular  dynamic 
simulation  code  based  on  Thee’s  MD  program.  Therefore  the 
initial  efforts  to  ascertain  the  cause  of  the  shrinkage  were 
focused  on  confirming  the  problem;  that  is,  insuring  that 
the  ca 1 cu 1  at i ona 1  cell  was  indeed  trying  to  shrink  and  that 
the  observed  behavior  was  not  due  to  a  problem  in  the 
computer  code.  If  the  shrinkage  were  caused  by  an  error  in 
the  ( N , P , H )  simulation  portion  of  the  code,  repeating  the 
principal  runs  using  (N,V,E)  simulation  should  show  a 
positive  pressure,  indicating  the  system’s  desire  to  expand. 
On  the  other  hand,  if  the  <N,V,E)  simulation  showed  a 
negative  pressure,  the  inclination  toward  smaller  volume 
would  be  confirmed. 


j 

S 


» 
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4  1 


To  test  this  hypothesis,  a  series  of  constant  volume 
<N,V,E)  runs  were  completed  at  0,  300,  and  1200*  K  using  a 
range  of  lattice  constants.  The  runs  consisted  of  1000 

scaling  steps  for  the  300°  runs  and  3000  scaling  steps  for 

_4 

the  1200*  runs  (2  x  10  psec  time  step).  The  initial  atom 
positions  and  velocities  for  the  300*  runs  were  the  final 
values  from  the  300*  principal  run  (Run  1  >  ,  but  the  1200* 
runs  and  the  0°  runs  were  begun  from  normal  lattice 
positions  with  appropr iately  scaled  random  velocities. 
Temperature  scaling  was  applied  each  time  step  to  insure  the 
final  temperatures  would  be  as  close  as  possible  to  the 
desired  value  so  that  the  total  energies  could  be  compared. 

The  results  of  these  runs  are  shown  in  Figures  15  and 

L 

f  16.  The  comparison  of  the  final  average  pressures  of  the 

runs  confirm  the  shrinkage  displayed  by  the  (N,P,H> 
principal  runs.  The  final  average  total  energies  of  many  of 
the  runs  were  too  similar  to  allow  conclusions  to  be  drawn, 
even  with  the  statistical  base  from  3000  time  steps.  It  was 
determined  that  statistically  significant  energy  differences 
would  require  unacceptably  long  runs,  so  no  further  attempts 
were  made  to  determine  the  resultant  energy  differences. 


l9§f  1°0  i.0fo  of  Shrinkage 

In  the  Tersoff  potential,  the  z  term  is  a  complicated 
combination  of  environmental  details  that  effect  the  bond 
strength.  Differences  in  bond  lengths,  bonding  angles,  and 
total  number  of  bonds  all  appear  in  the  calculation  of  z 
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Llncr^y  ( HJ  /k^ )  I  iO  I  Pressure  (GPa) 


.ire  15:  Final  Average  Pressure  vs.  Lattice  Constant 
for  0,  300,  and  1200  Degree  (N,V,E)  Runs 
<solid  symbols  annotate  N,P,H  result) 


Figure  16: 


Final  Average  Total  Energy  vs.  Lattice  Constart 
for  0,  300,  and  1200  Degree  (N,V,E)  Runs 
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In  an  effort  to  determine  the  cause  of  the  shrinkage 


phenomenon,  these  factors  were  examined  individually  in  a 
series  of  tests  that  will  be  discussed  in  this  section. 

Difference  fn  Bond  Lengths.  The  difference  in  the 
length  of  the  i-j  bond  compared  to  the  i-k  bond  appears 
cubed  in  an  exponent.  As  a  result,  one  would  expect  this 
term  to  play  a  major  role  in  determining  the  value  of  . 

To  test  the  importance  of  this  term,  a  (N,V,E)  simulation  at 
600  degrees  was  performed  with  X^  =  0,  which  effectively 
removes  the  difference  in  bond  lengths  from  the  z^j 
calculation.  The  run  consisted  of  500  scaling  steps,  100 

-4 

resting  steps,  and  4000  property  steps  with  2  x  10  psec 
time  step.  The  results  were  compared  with  an  identical  run 
made  with  X^  =  3.2394.  The  final  average  pressure  of  both 
runs  were  identical  (-0.59  GPa ) ,  and  the  final  average  total 
energies  showed  only  a  0.2*/.  difference.  The  small  changes 
noted  indicate  that  the  differences  in  bond  lengths  must 
remain  fairly  small  in  the  course  of  the  simulation,  so  that 
removal  of  the  term  in  an  exponent  has  negligible  effect. 

The  bond  length  difference  does  not  appear  to  be  the  driving 


force  for  the  observed  unrealistic  behavior. 


0D9yi§E  Degendence  of  Bond  Strength.  The  effect  of  the 
angle  defined  by  the  competing  bonds  was  examined  by 
calculating  the  value  of  z  for  a  single  j-i-k  triplet  and 
the  cor respond i ng  i-j  bond  energy  for  bonding  angles  between 
0  -  180  degrees.  The  result  is  a  surprisingly  flat  curve 

over  a  wide  range  of  bonding  angles  (Figure  17).  Caution 
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Figure  17;  Dependence  of  Single  Bond  Energy 
on  Bonding  Angle 


must  be  applied  when  estimating  the  importance  of  this 
result,  however.  Figure  17  shows  only  the  change  in  a 
single  bond  energy  as  a  result  of  a  single  bond  angle.  In  a 
tetr ahedr a  1 1 y  bonded  material  (silicon  in  a  diamond 
lattice),  change  in  any  single  bond  angle  effects  the  other 
bond  angles,  and  therefore  the  bond  energies.  The  result  of 
this  test  does  reinforce  earlier  statements  concerning 
Tersoff’s  dependence  on  coordination  number  rather  than  bond 
angle  to  stabilize  the  diamond  structure.  In  addition,  the 
results  also  hint  that  the  relative  insensitivity  to  bond 
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Figure  18:  Hypothetical  Two-Dimensional  Pair  Potential 
(arrows  show  direction  of  restoring  force) 


$ 


angle  may  be  part  of  the  problem.  This  possibility  will  be 
discussed  further  in  the  final  subsection  on  the  unexplored 
possib i 1 i t ies . 

Effect  of  Constraining  One  Bondiength.  One  final 
series  of  tests  was  completed  to  examine  the  effect  on  the 
total  energy  and  volume  of  the  ca 1 cu 1  at i ona 1  cell  when  one 
single  bond  length  within  the  cell  is  constrained  to  a  given 
value.  Before  discussing  the  results  of  this  test,  a  brief 
explanation  of  the  method  of  contraining  one  bondiength 
is  in  order . 

To  understand  the  atom  constraint,  consider  a  two- 
dimensional  parabolic  pair  potential  centered  on  some  rest 
separation  distance  between  two  atoms.  If  an  atom  is 
displaced  from  this  rest  position,  it  experiences  a  force 
attempting  to  return  it  to  the  rest  distance  (Figure  10). 
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19:  Hypothetical  Two-Dimensional  Pair  Potential  with 

Superimposed  Artificial  Energy  (no  restoring  force) 


If  an  artificial  parabolic  energy  curve  centered  on  the 


desired  position  is  added  to  the  pair  potential,  the  atom  no 


longer  experiences  a  restoring  force  and  the  bond  distance 


will  be  maintained  (Figure  19).  This  is  in  essence  th« 


technique  used  in  these  tests.  The  bond  length  between  two 


of  the  atoms  is  thus  constrained,  and  the  other  atoms  in  the 


computational  cell  are  allowed  to  "relax"  to  find  the  lowest 


energy  state.  The  artificial  energy  is  then  removed  and  the 


energy  is  calculated.  The  mechanics  for  this  atomic 


restraint  method  are  outlined  in  Sabochick’s  dissertation, 


and  interested  readers  are  referred  there  for  more 


information  on  the  technique  (14). 


This  series  of  tests  was  performed  for  a  range  of 


restrained  bond  lengths  and  a  variety  of  lattice  constants. 
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Table  V:  Minimum  Energy  <  eV )  of  a  64  Atom  System  with  a 
Single  Constrained  Bond 


Lattice  Constant  (Angstroms) 

Constrained 
Bond  Length 


( Angstroms ) 

5.41 

5.42 

5.4312 

5.44 

5.45 

2.117 

-296.007 

-296.037 

-296.041 

-296.024 

-295.981 

2.234 

-296.202 

-296.237 

-296.248 

-296.236 

-296.199 

2.352 

-296.250 

-296.291 

-296.307 

-296.298 

-296.266 

2.469 

-296.200 

-296.242 

-296.261 

-296.254 

-296.224 

2.587 

-296. 1 12 

-296. 159 

-296.183 

-296.180 

-296.155 

As  one  atom  in  the  ca 1 cu 1  at iona 1  cell  is  moved  closer  to  a 
neighboring  atom,  the  moved  atom  should  draw  its  other 
neighbors  with  it,  causing  a  slight  decrease  in  the  volume. 
Conversely,  as  two  atoms  are  held  slightly  apart,  a  slight 
increase  in  the  overall  volume  is  expected.  The  data  in 
Table  V  shows  the  model  behaving  properly,  with  a  slight 
decrease  in  the  lattice  constant  corresponding  to  the  lowest 
energy  for  constrained  bond  lengths  less  than  the  normal 
2.352  Angstroms  and  a  slight  increase  in  the  minimum  energy 
lattice  constant  for  constrained  bond  lengths  greater  than 
2.352  Angstroms.  This  tendency  is  more  apparent  in  Figure 
20,  where  a  spline  fit  to  the  data  for  the  constrained  bond 
lengths  2.117,  2.352,  and  2.587  Angstroms  shows  the  minimum 
energy  lattice  constants  to  be  approximately  5.427,  5.431, 
and  5.435  Angstroms,  respectively.  Therefore  the  shrinkage 
does  not  appear  to  be  a  result  of  an  erroneous  response  by 
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Figure  20:  Internal  Energy  vs.  Lattice  Constant  with 
Constraint  of  a  Single  Bond  Length 
(64  atom  system) 


the  secondary  neighbors  to  a  change  in  separation  distance 
between  two  bonded  atoms. 

yD§-£EioEed  Poss,ibi._lj.t_ies.  While  the  evaluation  tests 
discussed  in  this  chapter  failed  to  positively  identify  the 
cause  of  the  observed  shrinkage,  they  are  by  no  means 
exhaustive.  The  most  promising  clue  to  the  weakness  in  the 
potential  lies  in  the  relative  unresponsiveness  of  the  bond 
energy  to  the  bond  angle,  as  discussed  previously.  A  review 
of  the  role  of  the  z  term  will  help  explain  one  possible 
cause  of  the  shrinkage. 


term  will  help  explain  one  possible 


Table  VI:  Coordination  Number,  z. . ,  and  Minimum  Energy 
Separation  Distance  for^ Diamond ,  Simple  Cubic 
BCC ,  and  FCC  Structures 


Structure 


Coord  1 nat ion 
Number 


D 1 amond 
Simple  Cube 


3.4373 

2.3518 

6.0865 

2.5007 

8.9611 

2.6018 

14.628 

2.7298 

As  pointed  out  in  Chapter  II,  the  bond  strength  is 

modified  by  the  bonding  environment  through  a  changing  B 

value.  This  changing  B^  value  is  driven  by  a  z^j  value 

that  is  determined  by  the  surrounding  environment.  For 

periodic  structures  where  all  bond  lengths  are  equal,  it  is 

fairly  simple  to  calculate  the  resulting  z  value.  Table 

VI  shows  the  coordination  number  and  z  value  associated 

with  the  diamond,  simple  cubic,  BCC,  and  FCC  structures. 

By  taking  the  derivative  of  the  Tersoff  potential  with 

respect  to  the  bond  length  r^j  while  holding  z^j  constant, 

the  minimum  energy  bond  length  for  each  value  of  z . .  can  be 

J 

calculated  from 


These  minimum  bond  lengths  are  recorded  in  Table  VI.  From 


equation  (19)  and  the  known  relationship  between  z . .  and 


developed  in  Chapter  II,  it  is  apparent  that  a  smaller 


value  of  z..  equates  to  a  larger  B  and  a  therefore  smaller 

J-  J  -L  J 


minimum  energy  bond  distance.  Said  more  simply,  the  smaller 


z  becomes,  the  more  compact  the  system  becomes.  This 
-**  J 


offers  a  possible  explanation  for  the  observed  shrinkage,  as 
explained  below. 

Consider  a  system  containing  four  atoms.  If.  the 
absence  of  directionally  favored  bonding,  the  atoms  would  be 
expected  to  occupy  a  planar  arrangement  with  bond  angles  of 
120  degrees  and  all  bonds  of  equal  length.  This  structure 


gives  a  z  value  of  2.3139  and  a  corresponding  r  of 

lj  min 


3.5172  Angstroms.  If  the  structure  is  changed  so  that  all 


three  bonds  are  mutually  orthogonal,  the  value  of  z^  drops 


to  2.0  and  r  becomes  2.3132  Angstroms,  thereby  allowing 

i  J 


the  system  to  shrink. 

In  a  similar  manner,  it  may  be  that  the  tetrahedral 
structure  with  identical  bond  lengths  is  not  the  structure 


corresponding  to  the  lowest  value  of  z 


U 


Because  of  the 


extreme  sensitivity  of  the  z. .  term  to  coordination  number, 

J 


it  is  unlikely  that  any  structure  containing  more  than  four 
nearest  neighbors  is  energetically  favored,  but  a  different 
combination  of  bond  angles  and  bond  distances  for  the  four 


neighbors  of  an  atom  may  give  a  lower  z^j  value.  The 


existance  of  such  a  structure  could  easily  be  responsible 
for  the  shrinkage  noted. 
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VI.  Conclusions  and  Recommendations 


QQD£iy§i2D5 

The  directional  covalent  bonding  nature  of  silicon  is 
very  difficult  to  model  accurately.  Three-body  potentials 
such  as  those  proposed  by  Stil linger  and  Weber  and  Biswas 
and  Hamann  have  shown  some  successful  applications,  but  the 
strong  dependency  on  the  diamond  structure  limits  their 
usefulness.  The  approach  used  in  the  new  silicon  potential 
proposed  by  Tersoff,  explicit  modeling  of  the  bond 
environment,  seems  to  be  a  more  widely  applicable  model. 

The  purpose  of  this  research  effort  was  to  evaluate  this  new 
potential  in  a  molecular  dynamic  study  of  solid  silicon. 

The  Tersoff  potential  accurately  reproduced  the  correct 
bonding  lengths  and  bond  energies  for  the  diamond,  simple 
cubic,  and  FCC  structures  of  silicon  in  a  static  simulation 
(0°  K).  However,  as  the  system  was  heated  in  a  flexible 
boundary  simulation,  the  volume  decreased.  This  behavior  is 
not  observed  in  real  silicon.  Subsequent  constant  volume 
simulations  confirmed  that  decreased  volume  yielded  lower 
energy  of  the  model  when  the  simulation  temperature  was 
1 ncreased . 

Investigative  runs  were  completed  in  an  effort 
to  isolate  the  cause  of  the  problem.  The  effect  of  unequal 
bond  lengths,  varying  bond  angles,  and  changing  the  length 
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of  a  single  bond  in  a  64  atom  system  were  examined,  but 
none  provided  definitive  identification  of  the  source  of  the 
problem.  The  results  of  the  bond  angle  tests  suggested  the 
possibility  that,  while  Tersoff’s  potential  accurately 
reproduces  a  most  stable  configuration  having  four  nearest 
neighbors,  the  tetrahedral  diamond  structure  may  not  be  the 
lowest  energy  configuration. 

Because  of  the  unrealistic  behavior  under  heating,  the 
current  Tersoff  silicon  potential  is  unusable  for  dynamic 
simulation. 

?§£°!I?!D®Dd§iiQD5  • 
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The  generality  of  the  Tersoff  potential  vis-a-vis  the 
SS.W  three-body  potential  makes  continued  investigation  of 
the  Tersoff  style  of  potential  worthwhile.  Two  possible 
avenues  of  continued  study  are  readily  identified.  The 
first  is  to  experiment  with  the  values  of  the  10  variables 
in  the  Tersoff  potential  by  fitting  the  potential  to  dynamic 
properties  such  as  the  thermal  expansion  coefficient. 
Additional  parameters  may  need  to  be  added  to  accomplish 
this  fitting.  Indeed,  simply  increasing  the  cut-off  radius 
slightly,  say  to  3.6  Angstroms  (still  less  than  the  second- 
nearest  neighbor  distance  in  the  diamond  structure),  might 
help  prevent  the  shrinkage.  Such  a  move  would  strengthen 
the  angular  dependence  by  limiting  the  distance  one  neighbor 
could  move  before  interacting  with,  and  subsequently  raising 
the  coordination  number  of,  another  neighboring  atom. 


A  second  possible  extension  would  be  to  use  the 


potential  as  it  currently  exists  to  perform  a  series  of 
melting/annealing  runs.  The  atom  positions  could  be  sav’d 
after  each  run  in  an  attempt  to  identify  low-energy 
conf igurat ions  that  might  be  responsible  for  the  observed 
shrinkage.  Dr.  Tersoff  has  expressed  interest  in  pursuing 
this  action. 


Throughout  the  course  of  this  research,  dialogue  was 
maintained  with  the  author  of  this  potential.  Continuing 
this  interactive  relationship  offers  an  ideal  opportunity 
for  follow-on  researchers  to  play  a  significant  role  in  the 
development  of  a  tool  long  awaited  in  the  field  -  a  general 


BeB§Ddi2S  Derivation  of  Forces 


The  force  Fn  on  any  atom  n  is  simply  the  negative 
derivative  of  the  internal  energy  U  with  respect  to  the 
position  of  atom  n,  or 


d 

F=  -  -  U  (20) 

n  dr 

n 

Although  the  vector  notation  has  been  omitted*  equation  (20) 
actually  represents  three  scalar  equations,  one  for  each 
directional  component.  Since  the  form  of  each  component 
equation  is  identical,  only  one  component  is  shown  in  this 
appendix.  The  reader  is  expected  to  remember  that  parallel 
equations  apply  to  the  other  components. 

While  equation  (20)  appears  simple  enough,  applying  it 
to  Tersoff ’s  potential  can  become  confusing  because  of  the 
annotation  method  used.  Therefore,  before  the  derivatives 
are  discussed,  an  explanation  of  the  terminology  used  in  the 
potential  might  be  beneficial. 

The  basic  geometry  of  the  bond  triad  used  in  the 
calculation  of  the  bonding  energy  is  shown  in  Figure  21. 

The  central  premise  of  the  Tersoff  potential  is  that  the 
strength  of  the  bond  between  any  atoms  i  and  j  is  affected 
by  the  proximity  to  i  of  another  atom  k.  Therefore  for  each 
i-j-k  bonding  triad,  there  are  positional  derivatives  and 


55 


vW_ 


•  -  V* • 


Figure  21:  Geometry  of  Bonding  Atoms  in  Tersoff’s 
Potent i a  1  < 18 ) 

consequent  forces  acting  on  each  of  the  three  atoms. 

Due  to  the  lack  of  symmetry  caused  by  the  third  body  (k) 
term,  all  three  derivatives  must  be  solved  for  separately. 

The  annotation  used  in  the  Tersoff  potential  may  cause 
some  confusion  in  identifying  the  atoms.  The  "i-j-k" 
designation  of  the  potential  must  not  be  confused  with  the 
"n"  designation  of  equation  < 20 ) .  No  single  atom  in  the 
system  is  just  an  i,  j,  or  k  atom;  rather,  all  atoms  are 
all  three,  depending  on  which  bond  is  being  considered. 
Figure  22  shows  the  dual  annotation  applied  to  the  six 
possible  bonding  triads  for  any  three  given  atom  positions. 
For  convenience,  atoms  1,  2,  and  3  were  selected.  Across 
the  top  row  of  triads,  atom  1  is  atom  i,  atom  j,  and  then 
atom  k,  respectively.  Therefore  when  the  force  on  atom  1  is 


Figure  22:  Six  Bonding  Triads  Formed  by  Three-Atom 
System 

calculated,  it  must  be  the  sum  of  contributions  from  all  the 
triads  involved,  so  that 


Fn  1  k 

n , j , k  dr 


d 

- U 


d 

+  Z  - u 

i , n , k  dr , 


j 


d 

Z  —  U 

i  ,  j  ,  n  dr. 


(21  ) 


where  the  first  summation  is  over  all  triads  where  n  is  atom 
i,  the  second  sum  when  n  is  atom  j,  and  the  third  when  n  is 
atom  k  . 

In  the  actual  code,  each  triad  is  considered  once,  the 


dr  ,  dr  ,  and  dr,  derivatives  are  evaluated,  and  the 

i  J  * 


resulting  forces  are  added  to  the  previous  forces  for  the 
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appropriate  atom.  Before  each  new  time  step  begins,  the 
forces  are  r e 1 n i t i a  1 i zed  to  zero. 

The  form  of  the  Tersoff  potential  is 


u  <  r  .  .  )  =  f  ( r . , ) [ Aexp ( -A, r . . )  -  B, ,exp (-A^r.  .  ) ] 

ij  C  1J  1  lj  Ij  2  Ij 


(3) 


r  <  R  -  D 


r(r  -  R) 


f  (r)  =<.5  -  .  5s  in C 

c 


]  R  -  D  <  r  <  R  +  D  (4) 


2D 


r  >  R  +  D 


B  —  R  (1  -t-  h  7U  )-1/(2n) 
Bu  "  Bo(1  +  bzij) 


(5) 


2«  fc<rlk’  9<e>  ^«P<-tX3<rlk  -  Ty))3) 


(6) 


g (0)  =  1  +  c2/d2  -  c 2 / C d 2+  (cos  9  +  h)2] 


(7) 


Successive  application  of  the  chain  rule  can  be  used  to  find 


d  d 

- U  = - U 


dr 

—  U 


♦  i  u  d.:n 


dr 


dr 


ij 


dr 


dr  dr 

lk  i 


d  d 

- U  = - U 

dr  dr 

J  ij 


dr 

- ii 


dr 


d  dr 

- U  — 

dr  dr 

Jk  j 


d  d 

- U  = - U 


dr 

___lk 


dr 


dr 


ik 


dr 


k 


---  u  --ik 

dr  dr 

jk  k 
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These  three  equations  can  be  simplified  somewhat  by 
reverting  to  component  annotation  and  recalling  that 
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J 

In  terms 

of  the 

x  component 

,  it  is  apparent  that 

!r_ii 

-  .  trAl  . 

X 

.  x 
dri 

.  X 

'"J 

ru 

and  that 

similar 

equat ions 

apply  to  the  other 

derivatives, 

giving 

_:L 

< 

rik 

and 

dr 

— Ik 

dr 

= - JJt  = 

rx 

— i- 

drx 

j 

dr  x 
k 

r 

Jk 

Therefore  only  the  three  potential  derivatives  remain  to  be 
determined . 

Before  proceeding  with  the  derivatives  of  U  directly, 
it  might  be  more  convenient  to  get  rid  of  the  theta  term. 


ft 

; 


Using  the  Law  of  Cosines,  the  angle  can  be  related  to  the 
distances  separating  the  atoms  to  give 


cos  9  = 


2  2^2 

2rij  rik 


Taking  the  required  derivatives  gives 


d  cos  9 


•5C1/rik  ~  rik/rij 


d  cos  9 


=  .5Cl/r.  .  -  r.  /r?.  +  r?.  /(r.  .  r?.  )3 


ij  ik 


jk  lj  ik 


d  cos  9 

-  ~  rJk/<rij  rlk) 


Continuing  along  these  lines,  the  chain  rule  can  be  used  to 
determine  the  desired  derivatives  of  the  g  function,  giving 


dg  <  9) 


dg ( 9 )  d  cos  9 


d  cos  f 


with  similiar  functions  for  the  r  and  r  derivatives. 

ik  jk 

The  second  term  of  these  parallel  equations  is  already 
known,  and  the  first  term  is  simply 


?  5  22 

-  =  2c^ (cos  9  +  h)/Cd^  +  (cos  9  +  h)  ] 


>cvc 


Determining  the  derivative  with  respect  to  r 


is  simply  a 


matter  of  inserting  the  value  of  r into  the  equation  for  r 


and  solving. 


The  derivatives  of  are  somewhat  more  involved,  since 


B .  is  a  function  of  all  three  variables.  Once  again  the 
■*"  J 


chain  rule  is  applied,  yielding 


-  B 


d  d 

IT. 


drik  U 


Hr 


d  d  d 

_  q  =  — Q  _ — _  2 

drjk  1J  dzu  iJ  drjk  1J 


Again  the  first  term  in  all  three  equations  is  equivalent 


Solving  for  this  derivative  of  B  gives 

J 


d 

- B 

dz  ^  j 

ij 


B  (-<1  +  bzn  >“  l+l/(2n)  /2}  bz(n-D 
0  ij  ij 


Only  the  derivatives  of  z  remain  to  be  solved  for.  One 

ij 


final  time  the  chain  rule  is  applied  to  give 


2  ij  =  2 

dr.  .  J  k*i 


I  /f  (r 

i.i*  C 


)  expC-cA  <r  -  r  ) 3 - g(0) 

ij  3  ik  ij  dr 

ij 


+  f  (r  )  g<0)  exp{-cA  (r  -  r  )  33  > 
c  ik  3  ik  ij 


3Wrik-  ru  )j2} 


where  all  terms  on  the  right  side  of  the  equation  are 

included  in  the  summation  over  k.  Likewise  for  the  r  and 

ik 

r  derivatives, 
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i  *2  d 

2  lfc(rik:)  exP^~CAo(r..  -  r  >3J> - g<0> 

k*i,jl  C  11C  3  11C  1J  drik 

+  fc(rik)  9<®)  ®xp<-CA3(rlk-  r  >:3> 
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At  this  point  all  the  values  are  known,  and  back 
substitution  can  be  used  to  obtain  the  forces  on  each  atom 


of  a  triad 


When  these  forces  are  summed  over  all  the 


triads  an  atom  is  involved  in*  the  result  is  the  total  force 
experienced  by  that  atom. 
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Appendix  B:  Calculation  of  Virial  Term 


Kittel  defines  the  virial  term  used  in  the  pressure 
calculation  as  (10:223,224) 


VT  =  T  F  .  r 

i  i  i 


For  the  Tersoff  potential,  this  can  be  expressed  as 


VT  =  I  r  F1 

i  ,  j* i  ^  ^j k 
k*i  ,  j 


+  r  fJ 

J  ijk 


r  Fk 

k  ijk 


where  F~  is  the  force  on  atom  i  due  to  the  ijk  bond  triad 
1  J  K 


F-]  is  the  force  on  atom  j,  and  Fk  is  the  force  on  atom 

ijk  J  ijk 


k .  The  summation  is  then  over  all  bonding  triads.  Note 
that  vector  notation  has  been  repressed  for  convenience. 

This  may  look  like  a  more  complex  form,  but  it  can  be 
converted  to  a  form  that  is  independent  of  translation. 
This  is  accomplished  oy  application  of  Newton’s  third  law, 
according  to  which  the  total  forces  in  the  bonding  triad 
must  equal  zero,  or 


F^  +  F^  +  F^  =0 

ijk  ijk  ijk 


Multiplying  the  above  equation  by  r  ,  subtracting  the 


product  from  the  virial  term,  and  repeating  these  steps  for 
every  ijk  bonding  triad  gives 
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VT  =  Z 

i  ,  j*i 
k*i  >  j 


<  r  - 


r  )  F«J 
i  ijk 


(r  -  r  )  Fk 
k  i  ijk 


This  is  the  ca 1 cu 1  at i ona 1  form  used  in  the  MD  program. 
Notice  that  the  position  vectors  have  been  replaced  by  a 
difference,  which  removes  the  translation  dependence. 
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A  modified  silicon  potential  proposed  by  Tersoff  was 
evaluated  in  a  molecular  dynamic  simulation  of  diamond  silicon 
at  temperatures  well  below  the  melting  point.  Four  temperatures 
between  0  K  and  1200  K  wt  re  node  led  usin^  a  constant  enthalj.  y 
ensemble  of  CH  atoms. 

At  all  temperatures  above  0  K,  the  simulation  displayed 
increasing  atom  density  with  increasing  temperature.  This 
behavior  is  physically  unrealistic  and  makes  the  Tersoff 
potential  unsuitable  for  molecular  dynamic  simulation. 

A  variety  of  tests  were  performed  usinc  both  constan 
volume  and  constant  enthalpy  ensembles  in  an  effort  tc  de 
the  characteristic  of  the  Tersoff  potential  responsible  for 
thermally-induced  shrinkage,  but  results  were  negative. 
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